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Abstract 

This paper presents development of multi-input multi-output (MIMO) Generalized Pre- 
dictive Control (GPC) law and its application to reconfigurable control design in the event of 
actuator saturation. A Controlled Auto-Regressive Integrating Moving Average (CARIMA) 
model is used to describe the plant dynamics. The control law is derived using input-output 
description of the system and is also related to the state-space form of the model. The sta- 
bility of the GPC control law without reconfiguration is first established using Riccati-based 
approach and state-space formulation. A novel reconfiguration strategy is developed for the 
systems which have actuator redundancy and are faced with actuator saturation type fail- 
ure. An elegant reconfigurable control design is presented with stability proof. Several 
numerical examples are presented to demonstrate the application of various results. 


1 Introduction 

Over last decade, GPC has emerged as one of the leading control design strategies for robust 
control of dynamical systems. In early years (prior to late eighties) GPC’s applicability was 
essentially limited only to process control applications due to its demand on computational 
speeds. However, with the advances made in the computer technology over last decade 
computational speed is not a major concern for many real-life applications and control 
engineers have started using GPC for many main-stream applications. In recent years, 
GPC has become a viable alternative or in some cases even a preferred choice over well- 
known Hqo, H 2 , and /r-synthesis approaches. GPC has proved to be very effective when 
requirements on the robustness and performance are hard to achieve with traditional control 
designs. 

GPC belongs to a class of Model Predictive Control (MPC) methods. History of MPC 
dates back to late 7Q’s when the process industry showed keen interest in using these control 
methods. The control formulation at the time was mainly heuristic and algorithmic [1, 2], 
and exploited the increasing potential of digital processors. These controllers were closely 
related to~ the minimum timeoptimalc on t ro l me thodolo gy. The receding-horizonprincipln 
which is central to most of the MPC algorithms came about as early as 60’s [3]. As men- 
tioned earlier, MPCs became quite popular in the process industries where computational 
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speed was not a major concern. Also, many MPC algorithms were used on multivariable 
systems with constraints but no formal proofs of stability or robustness were available. 
Another parallel development took place using ideas from adaptive control which led to 
the development of self-tuning controllers [4] and extended horizon adaptive controllers 
(EH AC) [5]. This continued evolution of MPCs led to the emergence of the Generalized 
Predictive Control (GPC) methodology in late 80’s [16] which incorporates all major fea- 
tures of the predictive controllers in a unified framework. Various versions of the same 
common idea gave rise to the following different types of predictive controllers: Multi-step 
Multivariable Adaptive Control (MUSMAR)[6], Multipredictor Receding Horizon Adaptive 
Control (MURHAC) [7], Predictive Functional Control (PFC) [8], and Unified Predictive 
Control (UPC) [9]. 

MPC has also been formulated in the state-space setting [10], which not only allows 
the use of well established state-space theories for analysis but also provides the ease for 
extensions to multivariable systems. Moreover, it facilitates the use of stochastic theories 
and treatment of actuator/ sensor noise. The well developed estimation theory from state 
space methods can be easily incorporated without much complication. The perspective 
gained by working in these different domains made it possible to devise some simple tuning 
rules for ensuring stability and robustness for MPC systems. As a simple analogy, MPC 
controller can be viewed as an observer-based controller wherein its stability, performance, 
and robustness is determined by the observer dynamics, which can be fixed by adjustable 
parameters, and regulator dynamics, determined by MPC parameters such as weightings, 
horizon lengths, etc . In MPC, the control action at the current time instant is obtained by 
solring a finite horizon open-loop optimal control problem at each sampling instant. That 
is, the optimization problem yielding a control input is solved on-line at each time instant. 
The optimization process at each time instant yields an optimal control sequence and only 
the first control input in this sequence is applied to the plant. The novelty of MPC lies in 
its structure and not the control law itself. It essentially solves a standard nonlinear (or 
linear) optimal control problem. The fact that MPC solves the problem on-line and that it 
has the ability to naturally and explicitly handle the input and output constraints makes 
it a popular and powerful control paradigm. 

One major criticism of MPC methodology is that, except for some special cases, there 
is no general purpose theory which guarantees closed-loop stability of MPC. Although, in 
[11], some specific stability theorems are given for GPC using the state-space setting the 
general case stability results for GPC were lacking. Recently, in 90‘s, stability of GPC under 
end-point constraints was shown in [12, 13] where the equality constraints were imposed 
on the output after at the end of a finite horizon. A variation of this end-point constraint 
idea is used in this paper for establishing stability of Multi-Input Multi-Output (MIMO) 
reconfigurable control architecture using GPC. The stability proof uses Linear Quadratic 
Regulator (LQR) results and monotonicity properties of Riccati solutions. This paper also 
presents a clean and detailed formulation of GPC prediction equations and derivation for 
end-point constraint-based GPC control law for MIMO linear time-invariant (LTI) systems. 
A stable reconfiguration capability of the proposed GPC architecture is demonstrated using 
a simulation example. 
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2 GPC Control Paradigm 

GPC is the most generalized form of MPC. Among all different types of MPC schemes, GPC 
has maximum design freedom available for choosing design parameters. MPC design has 
several algorithms leading to different control schemes but all design schemes essentially use 
the same design paradigm. There are three basic elements of MPC: the prediction model, 
the cost function, and the control law. Every MPC uses some kind of plant model to predict 
future plant outputs over a pre-defined prediction horizon, N y = 7V 2 — Ni, where, N\ and 
N% represent lower and upper prediction horizons, respectively. The predicted outputs, 
y(t -f k\t) (k — Ni, . . . , JV 2 ) depend on the past inputs and outputs and on the future 
control signals u(t + k\t ), k = 1,2, ..., N u , where N u represents a control horizon. A 
set of future control signals is calculated by optimizing a suitable performance index so 
as to keep the plant output as close to the reference trajectory w(t + k) as possible. The 
performance index is typically a quadratic function of the predicted tracking error and 
control increments. The optimization is performed at each time step and the first element 
of the optimized control sequence is sent to the plant. The whole process is repeated again 
at the next time step after ‘"receding” the prediction horizon by one time-step. Because of 
this receding horizon feature of this control scheme it is sometimes referred to as Receding 
Horizon Control as well. The MPC control strategy is illustrated in Figure 1. Figure shows 
the predicted output based on plant model and predicted optimal control sequence based 
on predicted future error. In the figure, the control horizon N u is shown to be same as 
the prediction horizon N y , which is not necessary. In case of GPC, all design parameters, 
N \ , A 2 , JV U , and A ( control penalty in cost function) can be changed unlike other MPC 
schemes where one or more of these parameters are fixed. From this point on in the paper 
our focus will be on GPC control strategy unless otherwise mentioned. The control system 
block diagram for GPC is shown in Figure 2. 


3 GPC Control of SISO Systems 

-Consider- ay discrete-time-modeLd 
using a backward shift operator ( q ~ l ) as, 


A{q 1 )y(k) = q d B(q 1 )u(fc - 1) + C{q *) 


A 


SISO) system described 
(!) 
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Figure 2: MPC block diagram 

where u(k) and y(k) are the control and output sequences of the plant and £(k) is an 
uncorrelated random sequence. A , B , and C are the polynomials in the backward shift 
operator q -1 : 


A{q~ l ) 

— 1 T CLiq 1 + &2 q 2 * 

• • + a na q~ na 


= bo -+■ b\q ^ 4- 62*7 ^ 4- * 

• ’ + b n bq n>> 

Ciq- 1 ) 

= 1-1- ciq* 1 -f c 2 <?~ 2 4- • * 

■ + c nc q~ nc 


where d is the dead time of the system. A is the difference operator 1 - q~ l . This model 
is known as a Controlled Auto-Regressive Integrating Moving Average (CARIMA) model. 
For simplicity in the development, without loss of generality, C(q~ 1 ) is chosen to be 1. The 
plant (1) is then given by: 

A(q- 1 )y(k) = q- d B(q- 1 )u(k- 1) + ^ (2) 

The basic cost function used in GPC has the form 

n 2 n u 

J(NuN 2 , N u , A ) = £ [y(k + j) - w(k + j)] 2 + £ X(j)[Au(k + j- l)] 2 (3) 

j—Ni j = 1 

where, y(k + j) is an optimum j-step ahead prediction of the system output up to time 
k , w(k + j), j = 1 , 2 ,... is a future set-point or reference sequence. N± is the minimum 
prediction horizon, N 2 is the maximum prediction horizon, N u is the control horizon and 
A (k) is a control- weighting sequence. As seen in Eq. (3), the cost function is quadratic 
and penalizes future tracking errors over prediction horizon and control energy over control 
horizon. It is assumed that the control increments after control horizon are zero. The basic 
idea is to compute the optimal future control sequence, u*(k) = [u(k), u(k -f 1), .. .], such 
that the cost function J(Ni, N 2 , N u , A) is minimized. The control designer has to select 
the tuning parameters, N 2 , N u , and A, to meet certain stability and performance 
objectives. Once u*( k) is computed, only the first element of the sequence is used and the 
whole process is repeated at next time step. 

3.1 Prediction of future outputs 

In order to compute the cost function, which consists of future tracking errors, we need 
to compute the future outputs, y(k 4 - j) (N\ < j < N 2 ), using the best available plant 
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model (2). In case of SISO systems, these predictions can be easily computed using the 
Diophantine equation [9]. The development of prediction equations for SISO case using 
Diophantine approach is detailed in the Appendix 3.1. The prediction equations obtained 
in section 3.1 using Diophantine approach are very cumbersome to use and are not easily 
extendable to MIMO case. This motivates the development of state-space based formulation 
of prediction equations and control law. 


3.2 State-space formulation of prediction equations 

Consider a state-space description of the plant (Eq. 1) given as: 

x(Jfc + l) = Ax(ife) + BAv(Jfc) + B € f(*) (4) 

y(k) - Cx(fc) 4- £(k) 


The reason for using the state-space form of the system with A u(k) as input is that the 
algebraic complexity in the derivation of GPC control law reduces significantly in this form. 
The transformation of system equations with u(k) input to this form is given in detail in 
the appendix. In Eq. (4), the dimension of state vector is n — max(n a -fl, + d + 1, nc), 
and matrices A, B, Bf , and C are given by 


— ai 
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c d + 1 ~ Grf+i 


On &n 


C = [ 1 0 ••• 0 0 


One O-nc 


where a* are the coefficients of polynomial A , which is given by 

A(q ~ l ) — A A(q 1 ) = 1 + a\q 1 H — • 4- a n Q n + f- CLna+i<l (5) 

Note that, in matrix B, the d leading elements are zeros. The block dia^am representation 
of CARIMA model is shown in Fig. 3. If d = 0, the CARIMA model representation 
becomes that of Fig. 4. If in addition to d = 0, C(q~ 1 ) = 1, the CARIMA model takes the 
form of Fig. 5, For both cases of Figs. 4 and 5, we can get the corresponding state-space 
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Figure 3: The CARIMA model 
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A u(k) 



Figure 4: The CARIMA model with d — 0 


model. 

As mentioned previously, since the noise and disturbances in future are not known apriori 
for prediction of future outputs only deterministic part of the plant model is considered 
in the following development. The deterministic part of plant dynamics in polynomial 
description is given by 

A{q- l )&y{k) = B(q~ 1 )Au(k - 1) (6) 

Let the corresponding state-space model be 

x(k + 1) - Ax(k)+BAu(k) (7) 

y(k) - Cx(fc) 

Using the z-transform. we can obtain the z-domain transfer function H(z) as follows: 
H(z) = C(zl - A) -1 B 

c(i- 4r lB 

z 
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Figure 5: The CARIMA model with d = 0 and C(q 1 ) = 1 



CB CAB CA 2 B 
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By taking the inverse z-transform of the above equation, we get: 

y(k) = CBAu(fc - 1) + CABAu(k - 2) + CA 2 BAu(k - 3) + • • • 


The future outputs are then given by 

j/(fc + l) = CBAu(Jfc) + CABAu(fc - 1) 4- CA 2 BAu(fc - 2) + • • • 
y(k + 2) = CBAu(k + 1) + CABAu(k) + CA 2 BAu(k - 1) + ■ ■ ■ 

y(k + Ni) = CBAu(k + N 1 -l) + CABAu(k + Ni-2) + CA 2 BAu(k + Ni-3) + 


y{k + N 2 ) = CBAu(A + N 2 - 1) + CABAu(A; + N 2 - 2) + CA 2 BAu(A: + JV 2 - 3) + 
A general term y(k + j) (j = 1, 2, • • • , N 2 ) in the above equation can be written as 


y(k + j ) = ^CA^BA u(k+j-i) 

i—1 

j 


= ^ CA* -1 BAu(A; + j — i)+ CA i-1 BAu(fc + j- i) 

i = 1 i=j+l 


i=l 

J 


m = 0 
oo 


= CA i_1 BA u{k + j-i)+ ^2 CA m+j (x(fc - m) - Ax(fc - m - 1)) 


t=i 


m— 0 
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3 °° , N 

= Y CA l_1 BA u(k + j — i) + Y CA J (-A- m x(A: — m) - A m+1 x(fc — m — 1) J 

i=l 771=0 

j 

= Y CA i_1 BA u{k +j-i) + CA j x{k) 
i — 1 

Now using assumption of GPC control strategy that Au(fc+i) = 0 for i > N u , and denoting 


y(k + 1) 
y(k + 2) 


y = 


y(k + N!) 


y(k + N 2 ) 


and Au = 


Au(k) 
Au(k + 1) 


A u(k + N u - 1) 


we can write the predicted output as 


y 


QAu + T 


where 


CA 

CA 2 


T = 


CA Wl 


x(fc) 


CA W2 


and Q is given by 


CB 

CAB 


0 

CB 


0 

0 


CA Nl l B ca Wi 2 b CA Ni ~ 3 B 


CA JV2_1 B CA W2_2 B CA N2 ~ 3 B 


0 

0 

0 


CA N2 ~ Nu B 


( 8 ) 


Note that Eq. (8) allows easy computation of G compared to Diophantine-based technique. 
For SISO systems, it is interesting to note that the polynomials Ej and Fj in Diophantine 
equation solution have coefficients that relate to the A matrix of the state-space form of 
Eq. (4). To be specific, the elements in the first row of A? are the last N u coefficients of 
Ej + 1 in the decreasing order of the degree of q^ 1 and the first column has the coefficients 
of Fj in the increasing order of the degree of q~ l . 


3.3 Development of GPC for MIMO Systems 

In this section, we will develop GPC control law for the general case of MIMO LTI systems. 
For algebraic simplicity we will use state-space description for plant dynamics. Let the 
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state-space model for a general MIMO system with p inputs and q outputs is given by 

x(/c -hi) = Ax(fc) 4 - BAu(fc) 

y(fc) = Cx(fc) (9) 

where, x is the n-dimensional state vector, A, B and C are system matrices with dimensions 
n x n, n x p, and <7 x n, respectively. 

The cost function to be minimized is given by 

n 2 

J(Ni,N 2 ,N u ,X) = (y{k + j) - ™{k + j)) T {y(k + j) - w(k + j)) + 

j=Ni 

^ A(Au(fc + j — l)) T (Au(fc + j — 1)) (10) 

where, y is the g-vector of predicted output, w is the ^-vector of reference trajectory, and 
Au is the p ~ vector of input increments. In order to compute the optimal control input 
minimizing above performance function output predictions need to be computed over the 
prediction horizon. The following section develops prediction equations for MIMO case. 


3-4 Output prediction for MIMO Systems 

The procedure developed for predicting output in SISO case can be extended to MIMO 
systems as well with some added algebraic complexity. Consider a MIMO system of Eq. 
(9). We can compute the future outputs y(k 4 - j) ( j = 1 , 2 , • • • , N\, • • • , N 2 ) based on the 
plant information available up to time k using a recursion procedure. The predictions of 
state for future times starting at time k are given by 


x(fc 4 1) 
x(/c 4 2 ) 
x(fc 4 3) 


Ax(fc) 4 BAu(fe) 

Ax(A; 4 1) 4 BAu (k 4 1) = A 2 x(/c) ABAu(fc) 4 BAu(A; 4 1) 

A x(Jb 4 2) 4 BAu(A; 4 2 )= A 3 x(fc) 4 A 2 BAu(/c) 4 ABAu(fc 4 1) 
4 BAu(/c 4 2 ) 


x(fc 4 j) — A j x(k) 4 A - 7 1 1 BAu(/c 4 i) 

i—o 

The output at time k 4 j, y (k 4 j), is then given by 


y (k 4 j) = CA J’x(Jfc) 4 CAJ^BAu (k 4 i) 


Rewriting y(k 4 j) in the matrix form, we get 

( Au(fc) ^ 

Au(fc 4 1) 

~7(fc + j) = [ CA j l B CA J ^B • CAB CB ] r— 

Au(k + j — 2) 

\ Au(k+j - 1) / 

( 12 ) 
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Let us define 


and 


yi(k + Ni) 
y 2 (k + Ni) 


y = 


y(k + N 1 ) 
y (k + N 2 ) 


y q (k + Ni) 

yi(k + N 2 ) 
y 2 (k + N 2 ) 


Vq{k + N 2 ) 


(iV 2 -Ar 1+ i)g X i 


Aui(k) 

A u 2 (k ) 


Au = 


Au (k) 
Au (k 4- 1) 


Au (k 4- N u - 1) 


A u p (k) 
Aui(k 4- 1) 
Au2(k 4 * 1 ) 

A u p (k -f 1) 


A u\(k 4- N u — l) 
Av,2(k 4- N u — 1) 


A Up(k 4 - N u — l) 


J N u px 1 


(13) 


(14) 


Then predicted outputs are given in a compact form as 



CA Wl_1 B CA^“ 2 B CA Wi ~ 3 B 

0 


' CA Nl " 

y = 

CA W2_1 B CA N2 ~ 2 B CA W2 - 3 B 

. . CA N2 ~ Nu B 

Au 4- 

CA N2 


V ^ ' v ' 

G f 

y = GAu + f (15) 


Note that the form of G in MIMO case is same as that of SISO case (Eq. (8)) except the 
prediction horizon is from to The cost function of Eq. (10) can be rewritten in a 
compact form as 


J(Ni,N 2 ,N u ,\) = (y - w) r (y — w) + AAu r Au 

— . . - = (GAu + f — w)T(.GAu f.= w) i AAu T Au (16) 
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where, y, Au, G, and f are given by Eqs. (13)-(15), and w is given by 




wi(k + Ni) 

W2 (k + Ni) 

w (k + Ni) 
w (k + Ni + 1) 


w q {k + Ni ) 

w (fc + iV 2 ) 


wi(k + N 2 ) 
W2(k + N, 2 ) 



_ W q (k + N 2 ) _ 


Rewriting the cost function as, 


J = 


^Au t HAu + bAu 4- fo 


( 17 ) 


(18) 


where 


H = 2[G r G + AI] 
b = 2[(f-w) r G] 

^ = (f-w) T (f-w) 


and setting 


dJ 

SAu 


= 0 


we obtain the optimal control input A u* as follows: 


Au* = -H _1 b T = -(G r G + AI)- 1 [G r (f- w)] 


(19) 


The actual control signal sent to the process is the first p elements of the vector Au*. 


4 Stability of GPC control law 

One of the major drawbacks of GPC control law presented in Eq. (19) is that the closed- 
loop stability under this control law can not be guaranteed. However, this difficulty can 
be overcome if the performance function in Eq. (10) used to obtain the control law can be 
modified to include the penalty on the end-state of the system and GPC design parameters 
can be chosen in a certain way. 

The approach that is taken here to prove the stability of the GPC control law is to use 
the stability properties of the corresponding LQR problem. So, the optimal control problem 
can now be stated as follows: 

For th e LTI syst e m 


x(fc -hi) = Ax(fc) + BAu(fc) 
y(k) = Cx(fc) 


( 20 ) 
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determine the optimal control law that minimizes the cost with the end-point state weight- 
ing 

J(Ni,N 2 ,N u ,Q,X) = (x(fc + N 2 ) ~ w x (k + N 2 )) T Q(*(k + N 2 ) -w x (k + N 2 )) 

n 2 

+ E [y( fc + i) - w (^ + i)] T [y(^ + i) - w(fc + j)\ 

j=Ni 

N u 

+ E A(j)[Au(fe + j - 1 )] T [Au(/c Hh j - 1)] (21) 

3 = 1 

where, w x (fc + iV^) is the desired value of the state at the end of the prediction horizon, 
JVi = 1, N 2 = iV u = TV, Q > 0 and A(?) = A > 0. 

Since the stability properties of the system are independent of the reference input, we 
can ignore the reference input (or equivalently assume to be zero), i.e., set w x (/c + iV 2 ) = 0, 
and w (k + j) = 0, V?. Then, the cost function will be modified to 

J N = x T (k + N)Qx(k + N) 

N 

+ Y { y ( fc + j) T y( k +j) + ^ Au (* + j - i) Tu (k + j - i)j 

j = 1 

= x T {k + N)(Q + C T C)x(fc + N) 

AT- 1 

+ Y { x ( fc + j) r C r Cx(fc + j) + AAu (k + j) T Au(k 4- j) j 

3=0 

-x T (k)C T Cx(k) (22) 


4.1 GPC control law with end-point weighting 

The GPC control law derived in Section 3.3 is no more optimal for the cost function in 
Eq. (22) since this cost function now includes additional penalty on the end-state of the 
prediction horizon which was absent in the cost of Eq. (16). The GPC control law for this 
modified cost is derived next. 

The state evolution equation for the plant in Eq. (20) is given by 


x(fe 4- N) 


N - 1 

A N x(k) + Y A^ i_1 BAu (fc + i) 

i— 0 

/ 

A N x(k) + [ A n_1 B A w - 2 B • • • B ] 

^ A 


Au(fc) 
Au(A; + 1) 

Au(fc + N - 1) 

V 

AU 


\ 


/ 


= A N x(k) + CAu (23) 

Using Eqs! (23) and (15) the cost function in Eq. (22) can be re-written in more compact 
form as 


J N = iAu r [H + 2C r QC] Au + [b + 2x T (fc)(A JV ) T QC'] Au + f 0 + x T (fc)(A Af ) T QA N x(fc)(24) 
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where, 


H = 2[G T G + Al] 

b = 2f r G 
fo = ft 


Matrices G and f axe as defined in Eq. (15) and f is given as 

CA 1 
CA 2 


f = 


CA 


N 


x(fc) — Lx.(k) 


Let the cost function in Eq. (24) be further simplified to 

J N = iAu r HAu + bAu + f 0 (25) 

where 

H = 2 [G T G + AI + C r QC] 
b = 2x T (fc) [L T G + (A Ar ) T QC] 
f 0 = x r (fc) (L T L+(A N fQA N )x(k) 

Now the optimal control Au* can be obtained by setting the gradient of Eq. (25) to zero 
with respect to Au. That is, 

a a^ = 0 ^ 


Au* = -H" 1 ^ 

= - (G T G + AI + C r QC) _1 (G T L + C r QA w j x(fc) 

= — Kx(ft) (26) 


The following lemmas give the stability result for the corresponding LQR problem which 
will be used in proving the stability of GPC control law and the connection between GPC 
and LQR control laws. 


Lemma 1 For the system (20) the optimal control law minimizing the LQ performance 
index (22) is given by 

Au*(Jb) = -K {k)x{k) (27) 

where , gain K(k ) is given by 

K(k + j) = [B T P(k +j + 1 )B + XI\~ 1 B T P(k +j + l)A (28) 

and P is the solution of the following Riccati Difference Equation (RDE) 

T(k+j) = A P(fc+j+i)A-A T P(^+j+l)BrS r -^j+l)B+A^- 1 B r F(L+ J +l)A+^C 

P(k + N) = Q + C T C (30) 


with the boundary condition 
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Proof 1 Proof is given in the appendix. 

Lemma 2 For the system shown in Eq. (20) and the performance function given in Eq. 
(22), the GPC control law given by Eq. (26) is same as the LQR control law (Eq. C.26) 
evaluated at the same time instant 

Proof 2 Proof is given in the appendix. 

4.2 Stability result 

The RDE in Eq. (C.28) can be re-written by reversing the time index by defining 
Pm P(k + N — m), where m = 0, 1, 2, • • • , N — 1 
Then, the RDE in the forward form becomes 

P m + 1 = A T P m A - A T P m B(B T P m B + AI)- 1 B T P m A + C T C (31) 

with initial condition P 0 = Q + C T C. 

Now, if we consider the steady-state solution for the infinite horizon problem, i.e., let 
N — ► oo and let P m = P m +\ — P —constant as m — > oc, then, Eq. (31) becomes the 
Algebraic Riccati Equation (ARE) in variable P: 

P = A t PA - A T PB(B T PB + AI) _1 B t PA + C T C (32) 

The following lemma is state without proof which gives the conditions under which the 
matrix P has the stabilizing property. 

Lemma 3 If the system ( A , B , C) is stabilizable and detectable, and Po >z 0, then the 
matrix sequence, {Pm}™ = o> generated by the RDE (Eq. (31)) has a unique limit P which 
stabilizes the system and satisfies the Algebraic Riccati Equation (ARE). (Eq. (32)). 

Theorem 1 If the system given by Eq. (20) is stabilizable and detectable , Po h Pi, Q h 0? 
A > 0, N\ = 1, and N 2 — N u = N, then the GPC control law of Eq. (26) stabilizes the 
system. 

Proof 3 In order to show stability of GPC control law, it suffices to show that the feedback 
gain K used in Eq. (26) is stabilizing. From Lemma 2, this is equivalent to showing that 
the corresponding LQR gain in Eq. (C.26) is stabilizing. Note that as long as the GPC 
parameters (N±, N 2 , N u , and X) are fixed the GPC control is equivalent to a constant-gain 
state feedback. 

Rewrite the RDE (Eq. (31)) as a set of two equations: 

Pm ~ A 1 Pm A - A 1 P m B{B x P m B 4- A T)~ l B T P m A -f Q(rn) 

Q(m ): - C T C+ (P m - P m +i) (33) 

where, first of Eq. (33) is nothing but an ARE for a given m. They are also referred to 
as Fake Alg ebraic Riccati Equations (FARE) in the~ literature: Now,~by Lemma~3 if the 
system is stabilizable and detectable and Pm — P m +i 0 , Pm wM stabilize the system, i.e., 
~A — A — B(B T P m B+ X I)~ l B T P m A will have all of its eigenvalues strictly within the unit 
circle. Next, we will use the monotonicity properties of the RDE and the property of FARE 
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to show the stability of GPC. 

Consider matrix Pk defined as follows 

Pk = Pfc-i — Pk 

Then , from the result (Fact 5.2.1) on Page 95 of [22], Pk satisfies the RDE 

P k+ i = A T k P k A k - A T k P k B(B T P k B + X k )- 1 B r P k A k 


where 


Now, if 


A k = A-B(B T P k B+XI)- 1 B T P k A 
Q(k) = CFC- (fc= 0 
X k = B r P k B + A J >: 0 


Pi = P 0 - Pi h 0 


by the property of the solution of the RDE , 


P N = Pn-i ~ Pn bO. 

Then, from first of Eq. (S3), Pn - i stabilizes the system. Since , Pn- i = P(k -f 1), the 
LQR controller gain K(k) in Eq. (C.26) stabilizes the system. Now, from Lemma 2, GPC 
controller gain K in (Eq. (26)) stabilizes the system ■ 

4.3 GPC control law for non-zero reference trajectory 

In the case when the desired trajectory is non-zero time varying the GPC control law (Eq. 
(26)) takes slightly different form to account for the tracking error. Consider the LTI system 

x(fc +1) = Ax(fc) 4- BAu(fc) 
y(fc) = Cx(fc) 

and the cost function 

J(Ni , N 2 , N u , A) = (x(fc 4- N 2 ) - w x {k + N 2 )) T Q(x(k 4- N 2 ) - w x (fc 4- N 2 )) 

n 2 

+ E [y(* + j) - w (* + j)] t [ y( k + j) - w ( fc + j)] 

j=N x 

+ E A0)[ Au (fc + j - l)] r [Au(fc + j — 1)] (34) 

i = i 

where w x (A; 4- N 2 )) is the desired value of x at t = k + N 2 . 

J = (CAu + A N2 x(fc) - w x (k + 7V 2 ))) r Q(CAu 4- A N *x{k) - w x (k + N 2 ))) 

+(y - w) T (y - w) + AAu r Au 
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Substituting y by from Eq. (13), we can rewrite the cost function as 

J = hu r HAu + bAu + f 0 

z 

where 

H = 2(G t G + AI + CQG) 

b = 2[(A N *x{k)-w x (k + N 2 )) T QC + (f-w) T G\ (35) 

f 0 = (A^ 2 x(fc) - w x (k 4- N 2 )) T Q(A N2 x(k) - w x (k + N 2 )) + (f- w) T (f- w) 

The optimal control law is then given by 

Au* = — H -1 b T . (36) 

where H and b are given by Eq. (35). 

4.4 Computation of desired state trajectories 

In most real life situations it is safe to assume that the desired output trajectory is given 
to the control designer. However, the control law computation (Eq. 36) requires that the 
desired state trajectories are known. This section gives the procedure for computing the 
desired state trajectories given the desired output trajectories for a general case. 

Consider MIMO LTI system given by the following state-space description 

x(fc + 1) = Ax(fc) 4 BAu(fc). (37) 

Let w (k 4* AT 2 ) denote the desired output trajectory that needs to be tracked by the sys- 
tem and let w x (k 4- A 2 ) denote the corresponding trajectories for the desired state. Then, 
w x (k 4 7V 2 ) must satisfy the following conditions: 

w (k + N 2 ) - Cw X (k + N 2 ) (38) 

(A — T)v/ X (k + JV 2 ) = 0 (39) 

The first condition (Eq. 38) is a direct consequence of output equation. The second condi- 
tion (Eq. 39) is derived next. 

One consideration in tracking problem is that the choice of w x (k + iV 2 ) should be con- 
sistent with the zero steady-state output tracking error. Consider the GPC control law 
with tracking error feedback: 

Au*(A;) = — Ke(fc) (40) 

where 

e(k) = x(/s) — w x (k) 

The closed-loop error dynamics for the system given by Eq. (37) with control law (40) is 
given by 

e(k 4- 1) = x(fc 4- 1) - w x (ft + 1) 

= Ax(fc) + BAu* (A:) - w x (k + 1) 

— (A - BK)e(fc) 4 Aw x (k) - w x (A; 4 1). 
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At steady-state, for zero tracking error we need 

e(k 4- 1) = e(k) = 0, and w x (fc) = w x (h 4 1). 
Then, to satisfy the end-point constraint of zero tracking error, we need 


(A - I)w x (k 4 7V 2 ) = 0. 
Combining the conditions of Eqs. (38) and (39), we get 


A -I 
C 


w x (fc + N 2 ) 


0 

w (k + N 2 ) 


( 41 ) 


Note that Eq. (41) typically represents under-determined system of equations as matrix on 
the left hand side is not generally invertible. In such as case, a least square solution can be 
obtained for the desired state trajectories as 


w x (/c 4 -N 2 ) 


( 

A — I 

v 

o 


C 

) 

[ w (k + N 2 ) 


(42) 


where (•)! denotes pseudo-inverse of (*). 


4.4.1 Special Case: System with redundant actuators 

For the system which has more inputs than outputs, i.e. p > q, the transformation of the 
system with u (k) as input to the form which has Au(fc) as input, the detectability may be 
lost. To see this, consider the matrix 


AI-A 
C 

where A and C are defined as in Eq. (B.16). When A = 1, rank of (AI — A) is less than 
n — p. The rank of C can at the most be q. So, the matrix in Eq. (43) must have rank less 
than n, which is not a full column rank. That is, the detectability is lost. However, this 
problem can be addressed as follows: 

Consider a system with Au (k) as input and let the number of inputs be greater than 
the number of outputs, i.e., p > q. 


x(fc + 1) = Ax(fc) 4 B Au(&) 

y(fc) = Cx(fc) (44) 


Now consider a modified output y(k) such that 





y(k) = Cx(A;) 

(45) 

where 





y(*) = 

y(*0 

u(* - 1) 

, c = 

l 1 

,P 

X 

A o 
•r 

X 

i i 

(46) 
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It can be verified that with respect to this new output the system is detectable. Thus, the 
modified system now takes the form: 

x(/c -f 1) = A x(/c) -f BAu (k) (47) 

y(k) = Cx(fc) (48) 

which is stabilizable and detectable. For this new system, let the cost function be modified 
as follows: 


N 


N 


J = x T (k + N)Qx(k + N) + Y^ y T ( k + j)Qy( k + j) + Yl + j ~ 1 ) 


(49) 


3 = 1 


3=1 


where 



(50) 


Now, it can be seen that minimizing the cost function given by Eq. (49) for the system 
given by Eq. (47) is same as minimizing the original cost function (Eq. 22) for original 
system (Eq. 44). The advantage in dealing with transformed system is that it satisfies the 
conditions of stability theorem and therefore allows computation of GPC control law. 

Now, to show that the stability proof still holds consider the following: 

We have already shown that the modified cost function is essentially same as the original 
cost function. Now it remains to show that the Riccati equations used in Lemma 1 with 
appropriate boundary condition do not change as a result of transformed system. The 
equations in Lemma 1 for the modified system can be rewritten as: 


Au *(k) = “K (fc)x(fc) 

where 

K(k + j) = [B T P(k +j+ 1)B + AI]- 1 B t P(£; +j+ 1)A 
P(k + j ) = A T P(k + j + 1)[A — BK(k + j)} + C T QC 

with 

P{k + N ) - Q + C T QC. 

By direct substitution, it can be easily verified that 

C T QC = C T C. 

This essentially proves that the GPC control law derived for the transformed system stabi- 
lizes the original- system.- ~~ - - - 
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5 Reconfigurable control 

In this section, we will extend the stability result obtained in previous section to the case of 
reconfigurable control architecture. We present a reconfigurable control scheme wherein the 
system has some redundancy in control actuators and in the event of actuator failures, like 
saturation for example, the system can reconfigure itself to re-allocate the control signal to 
other set of actuators. 

Consider a strictly proper stabilizable and detectable MIMO LTI system from Eq. (9) 
with p inputs and q outputs. 


x(*+l) = Ax(fc) 4- BAu(fc) (51) 

y(k) = Cx(fc) 

Let us suppose that the system has some redundancy in actuators, i. e., p> q and some of 
the p inputs can control more than one outputs. The reconfiguration problem is defined as 
follows. If one or more of the p actuators saturate (i.e., become “inactive”) how to re-allocate 
the control signal to actuators which are still “active” such that the lost degrees of freedom 
due to saturation in controlling certain outputs are regained via re-allocation. Moreover, 
to ensure that the overall stability of the system is maintained under such reconfiguration. 
Before reconfiguration scheme is presented and stability of such scheme is established some 
definitions are in order. 


5.1 Definitions 

Definition 1 Reconfiguration matrix: Any q x p matrix with binary entries (0 or 1) is 
called as Reconfiguration Matrix , Q rc , and is used to set the control priorities for each 
actuator. 


Q rc is used to select which actuator is used to control a certain output. The size of the Q rc 
matrix is q (number of outputs) by p (number of inputs). Each element of Q rc is either 1 
or 0 indicating weather a particular actuator is allowed to control a specific output or not. 
For example, let us take a 5-input 2-output system. An example Q rc matrix could be 


Q 


rc — 


0 10 0 0 
0 0 110 


(52) 


In this example of Q rc , the first output can be controlled only by the second input, and the 
second output can be controlled by both the third and fourth inputs. 

This matrix is then used to “reconfigure” the H matrix (where H is the transfer function 
matrix of system in Eq. (51)) as follows which causes the reconfiguration. Reconfigured H 
matrix is denoted by H rc and is given by 


H rc = Q r c ® H (53) 

where, ® indicates the outer product. This outer product effectively modifies the B and 
The c matrix of the system to"account “for reconfiguratidn“ — — 

Definition 2 Let Q be the set of all allowable Q rc s such that the outer product in Eq. (53) 
produces a set of reconfigured stabilizable and detectable systems. 
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That is, 


Q = {Q* c |(A,B i ,C i) is stabilizable and detectable Vxj (54) 

where, B; is the reconfigured B matrix and C* is the reconfigured C matrix as a result of 
outer product of Q l rc with G. 

Definition 3 The matrix designed to prioritize the order of the actuators in which they 
are used in the event of saturation failure is called as Reconfiguration Order Matrix. 


The design of the reconfiguration matrix Q rc is based on the criteria designed suitably based 
on the application, for example one of the criteria could be the saturation of the actuators. 
The order in which the redundant actuators are used is based on the reconfiguration order 
matrix, Q ro . For example, for the Q rc defined in Eq. (52) a reconfiguration order matrix 
could be 


Qro 


4 12 2 3 
4 3 112 


(55) 


This Qro matrix informs the reconfiguration algorithm to use the second input to control 
the first output until it saturates followed by actuators 3 and 4 until they saturate followed 
by actuators 5 and then 1. Similar logic is used for the second output. Each time the 
control law is evaluated, the first actuators are used to calculate the control input. The 
input is tested for saturation, if it is saturated, the Q rc matrix is set to the next actuator 
(or set of actuators). This process is repeated until there is no change in the Q rc matrix or 
there are no more free actuators left. Then, the resulting control is sent to the plant. This 
process is repeated at each cycle. The repeating of this process ensures that the reverse 
order of actuators is used as the actuators unsaturate. 


5.2 Stability under reconfigurable control 

The stability of the reconfigurable control architecture described above will be established 
in two steps. In the first step, it will be shown that the set of plants obtained after 
reconfiguration is stabilized by the GPC control law of Eq. (26). Then, in the second 
step, it will be shown that the system remains stable during the switching of the control 
law under reconfiguration. The stability of the reconfigured system is established by the 
following theorem. 

Theorem 2 All systems resulting from reconfiguration of system given by Eq. (51) due to 
reconfiguration matrix Q rc are stabilized by the control law of Eq. (26) iff Q rc belongs to 
the set Q. 

Proof 4 Proof of this theorem is the direct consequence of Definition 2 and Theorem 1. 

Note that the control law of Eq. (26) stabilizes the original system (Eq. (51)) as well 
according to Theorem 1. Having established that both systems - system before reconfigu- 
ration and the system after reconfiguration - are stabilized by the control law of Eq. (26) , 
to complete the stability argument we need to show that the stability is retained during the 
process of reconfiguration. For that, let the system after reconfiguration is described by 

x(fc -hi) = Ax(fc) + BAu(fc) (56) 

y(k) = Cx(A;) 
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Let us suppose that the system (51) is under closed-loop control with control law given 
by Eq. (26). Suppose the system has actuator saturation at time instant k — r and is 
reconfigured by matrix Q rc € Q to yield the system of Eq. (56). That means, the dynamics 
of the system (51) starting at say initial time k = 0 is governed by the dynamics of system 
(56) after reconfiguration (i.e., for V k > r). Since the GPC control of Eq. (26) is stabilizing 
for system (51), x(r) is bounded. Let x(r) be the initial condition for the system (56) at 
time instant k = r. Since x(r) is related to x(r) by some similarity transformation, x(r) 
is also bounded. Now, for time k > r, system evolves according the dynamics of Eq. (56). 
Prom Theorem 2, the system (56) is also stabilized by the control law of Eq. (26), i.e., x(fc) 
remains bounded for all r < k < oc. Thus we have shown that the system remains stable 
for all 0 < k < oo under control law of Eq. (26) ■ 


6 Numerical Example 

In this section, a numerical example is given to demonstrate stable reconfiguration capability 
of the GPC control methodology presented in Section 5 using a short-period approximation 
model of an aircraft. Two different cases of reconfiguration are presented where each case 
can be considered to represent a different type of actuator limitation and/or failure. 

The short-period approximation model of an aircraft considered is: 


x(t) = 

0 -1.3677 ’ 

x(t) + 

' -0.0665 

-0.0029 

-0.0029 " 

1.0000 -1.5087 

-0.0862 

-0.0043 

-0.0043 

y(t) = | 

0 0.2500 ] x(t) 






where the inputs are elevator deflection, left aileron deflection, and right aileron deflection. 
The eigenvalues of the open-loop system are —0.7543 ± 0.8937L The maneuver under 
consideration is the step change in the pitch-rate. The initial actuator configuration is such 
that only the desirable input is the elevator input to accomplish this task. The left and 
right ailerons are considered as redundant actuators and are to be used if elevator fails. 
For demonstrating reconfiguration capability two different cases are considered: Case (i) 
system reaches the steady-state condition and elevator fails to hold its position and elevator 
input suddenly drops to zero; and Case (ii) before system reaches to target output elevator 
freezes in its position and some constant input is acting on the system but not adequate to 
reach the desired position. These two cases are chosen as they represent two quite different 
dynamic characteristics. In the first case the redundant actuators (both ailerons) have to 
compensate for no elevator input and try to return system to desired steady-state condition 
whereas in the second case the ailerons have to make-up for the inadequate control input 
to reach to the desired steady-state. Given below are reconfiguration results based on the 
methodology presented in Section 5. 

Reconfiguration strategy and GPC control law design: For both cases, the steps 
involved in the control law design are same and are given below. The reconfiguration matrix 
Q TC has its in itial settin gs as: 

Qrc = [ 1 0 0 ] 5 
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With this newly defined C matrix, the system is stabilizable and detectable. Now, choosing 
the GPC design parameters as: Nj = 1, iV 2 = N u = 5, A = 0.1, and end-point weighting 
matrix Q = diag{ 150 800 1} we can obtain a stabilizing control law. The Riccati solutions 
P 0 and Pi for these set of design parameters are given below. 

P 0 = Q + C T C 

' 150.0002 0.0009 0 

= 0.0009 800.0044 0 

0 0 1.0000 

' 147.0957 7.5878 0.1832 ' 

Pi = 7.5878 608.4186 0.7829 

0.1832 0.7829 0.0921 

These Riccati solutions satisfy the conditions of Theorem 1 as the eigenvalues of the differ- 
ence P 0 - Pi are 

Eig(P 0 - Pi) = [ 191.8935 2.6266 0.8780 ] > 0 P 0 -Pih 0. 

When the system is reconfigured (in cases (i) or (ii)), the reconfiguration matrix Q rc is 
given by: 

Qrc = [ 0 1 1 ] . (61) 
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Then, from Eq. (53) the system transfer matrix after reconfiguration effectively becomes 


*(*) = 

y(t) = 


0.0000 -1.3677 
1.0000 -1.5087 


x(i) + 
[ -0.0000 0.0313 ] x(t) 


— 0.0234 -0.0234 
-0.0345 -0.0345 


u(t) 


(62) 


Again, discretization of the system with sampling time, T — 0.05s, yields 


x(k + 1) 


0.9983 -0.0658 
0.0481 0.9257 


x(k) + 


- 0.0011 - 0.0011 
-0.0017 -0.0017 


u(k) 


y (*0 


-0.0000 0.0313 ] x(k) 


(63) 


Again, as before, using the transformation given in the Appendix B one can transform the 
above system to a new system with A u(k) as input as follows: 



' 0.9983 

-0.0658 

-0.0011 -0.0011 ' 


' -0.0011 

-0.0011 ' 

x(A; + l) = 

0.0481 

0.9257 

-0.0017 -0.0017 

x(k) + 

-0.0017 

-0.0017 

0 

0 

1.0000 0 

1.0000 

0 


0 

0 

0 1.0000 


0 

1.0000 

y(*0 = [ 

-0.0000 

0.0313 

o 

o 





A u(k) 


(64) 


This system again is not detectable and one can define new fictitious output as before and 
define a new C matrix as 


C 


aug ~ 


-0.0000 0.0313 0 0 

0 0 1.0000 0 
0 0 0 1.0000 


(65) 


The modified system with this new C matrix is again controllable and observable. Now if 
we define the new GPC parameters to be Nj = 1, N 2 — N u — 5, A = 0.1, and the end-point 
weighting matrix as Q = diag{ 130 200 1 1}, we obtain the following Riccati solutions 
which satisfy the stability condition of Theorem 1 as shown below: 


Po 


Pi 


Eig(P 0 — Pij = [ 


= Q 4- C T C 


130.0000 

-0.0000 

0 

0 


-0.0000 

200.0010 

0 

0 


0 

0 

1.0000 

0 


0 

0 

0 

1.0000 


129.9835 

0.2793 

-0.0146 

-0.0146 

0.2793 

171.7884 

-0.0275 

-0.0275 

-0.0146 

-0.0275 

0.0909 

0.0000 

-0.0146 

-0.0275 

0.0000 

0.0909 


28.2154 0.0133 0.9095 0.9091' > 0 4- Pq-P\ h 0. 


Thus the conditions in the Theorem 1 are satisfied and the GPC control law (Eq. (26)) 
stabilizes the system. 
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Case (i): Consider the first of the two cases mentioned above. In this case, it is assumed 
that after system output (pitch rate) reaches its steady-state value of 1.15 deg/sec the 
elevator input abruptly goes to zero (at t = 2 secs) (potentially this could be a result of 
possible failure) and system is reconfigured using Eq. (6). 

Figure 7 shows the system response before and after reconfiguration. It can be seen 
that the system reaches a zero steady-state error within 2 seconds and when elevator input 
abruptly goes to zero at t = 2 secs, the output starts drifting from the desired value. At 
t = 2 sec reconfiguration takes place and ailerons take over and bring back the output back 
to the desired value rendering zero steady-state error. The input time histories are given in 
Fig. 6. As shown in the figure, the elevator deflection has abruptly dropped to zero after 
reaching an approximate steady-state value of —1.8 deg and at the same instant ailerons 
take over and bring the output back to the desired value. The ailerons reach the steady 
deflection of about -19 deg. The cost and tracking error plots are given in Figs. 8 and 9, 
respectively. 



Figure 6: The GPC control law 



Figure 7: The Plant Output 
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Figure 9: The Plant Output Error 


Case (ii): In the second case, it is assumed that the elevator input freezes to some 
nonzero constant —6 deg before the system output (pitch rate) reaches its desired value 
(5.73 deg/sec). The simulation results for this case are shown in Figures 10-13. As shown 
in Fig. 10, the elevator input freezes approximately at —6 deg when the system output has 
reached only 90% of its final value. At this time ailerons take over as a result of reconfig- 
uration and ailerons deflection reaches a steady-state value of approximately —21 deg as 
pitch rate reaches to desired value of 5.73 deg/sec. As seen in Fig. 11, overall tracking 
performance is very satisfactory without steady-state error. The performance function and 
tracking error time history is shown in Figs. 12 and 13, respectively. 

7 — Co ncludi ng Remarks 

This paper presented a comprehensive development of the basic GPC control architecture 
and derivation of control law for SISO as well as MIMO systems. For SISO systems, both 
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Figure 10: The GPC control law 



Figure 11: The Plant Output 

z-domain (i.e., polynomial-based) as well as state-space based control law development 
was given. Stability of GPC based on end-point constraints was presented as a basis for 
establishing stability of a novel reconfigurable GPC control architecture. The stability 
conditions for reconfiguration scheme were presented. A numerical example consisting 
of longitudinal control of aircraft was presented to demonstrate the effectiveness of the 
proposed reconfiguration methodology. 
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Figure 12: The Cost Function 



Figure 13: The Plant Output Error 
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Appendix 

A Diophantine equation 

The Diophantine equation offers a convenient way to obtain a division of two polynomials 
and has the following form: 

1 = Ej(q~ 1 )A(q~ 1 ) + q^F^q -1 ) (A.l) 

where A(q~ l ) = AA(q~ 1 ). The degrees of polynomials Ej and Fj axe j — 1 and n a , 
respectively. This equation will be useful in writing predictions of future outputs. The 
Diophantine equation has a unique solution, Ej and Fj , given A(q~ 1 ) and the prediction 
interval j. 

Consider Eq. (2) and multiply it by EjAqi to obtain 

E j AA(q~ 1 )y(k + j ) = EjB{q~ l )Au{k +j-d- 1) 4- E£(k + 3 ) (A.2) 

Now substituting for EjAA from Eq. (A.l) gives 

y(k +j) = Fj(q~ 1 )y(k) + Ejq'^Biq-^Auik +j-d- 1) + Ejiq-'Mk + j) (A.3) 

As the degree of polynomial Ej(q~ 1 ) = j — 1 the noise terms in Eq. (A.3) are all in the 

future and therefore can be omitted from the prediction as any better prediction of noise 
is not possible. The best prediction of y(k -f j) is therefore given by 

y(k 4- j) = Gj{q~ 1 )Au(k + j-d- 1) + Fj{q- l )y{k) (A.4) 

where, G J (q~ 1 ) = Ej(q~ 1 )B(q~ l ). Since the degrees of Ej and Fj are j — 1 and n a , the jth 
prediction of output involves the past outputs: ( y(k),y(k — 1), * * • ,y(k — n a )) and inputs 
(A u(k + j — d - 1), Au(k H- j — d — 2), • • • , A u(k — d - rib)). The polynomials Ej and Fj 
obtained from Eq. (A.l) are given by 

Ejq- 1 ) = ejfl + ej,i q~ x + ■ ■ ■ + ejj-i 

Fj(q~ l ) = fj,o + fj, lQ 1 4 1- fj,naQ ™ 

Now, the same Diophantine equation can be used to obtain Ej+ \ and Fj+ 1 , where 

Fj+lfa" 1 ) = fj+l ,0 + fj+l,lQ 1 H b fj+l t naQ 

and the polynomial Ej+ \ is be given by 

Ej+\{q~ l ) = Ej(q~ 1 ) + e J+ i jq ~ 3 

with ej + ij = fjfi. The polynomial G J+ \ can now be obtained recursively as follows 

Gj + 1 = Ej+iB — {Ej + fjoq 3 )B 

= Gj + f jfi q- j B (A. 5) 
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where 


G i(q *) = 9j,o + 9j,iq 1 + 9j,2q 2 + • • ■ ( A -6) 

From Eqs. (A. 5) and (A.6), it can be seen that the first j coefficients of G j+i will be 
identical to those of Gj and the remaining ones are given by 

9j+ij+i = 9jj+i + fj.obi i = 0, • * * , nb 

In the case when system has a dead time of <2, the output of the system at any time instant 
k is the result of input applied to the system at time k — d. As a result, the horizons Ni 
and N 2 should be chosen such that N 2 > Ni > d 4- 1. 

Now consider the following set of j-step ahead optimal predictions for j — <2-f 1, * - • , iV 2 : 


y(k + d + 1) = G d+ iAu{k) + F d+1 y(k) 
y(k 4- d + 2) = G d j r 2 ^.u(k + 1) 4- F d+ 2 y(k) 


y(k + Ni) = G Nl Au{k + Ni -d- 1) + Fs 1 y(k) 

y{k + N 2 ) = G N2 Au(k + N 2 -d-l) + Fx 2 y(k) 
which can be re-written as: 

y = G Au + F {q~ l )y{k) + g\q~ l )Au{k - 1) 


(A- 7) 


where 


' y(k + d+ 1) ' 
y{k 4- d + 2) 

y(k + iVi) 

. y(k + N 2 ) . 


Au = 


Au(k) 

A u(k 4- 1) 


[ Au(k + N U — 1) 


go 0 0 - - • 0 

9i 9o 0 0 

9Ni-d-l 9Ni-d-2 9Ni-d-3 "• 0 


L 9N 2 -d- 1 9N 2 -d- 2 9N 2 -d- 3 • • • 9N 2 -N u -d J 



{G d+l {q l ) - g 0 )q 



( G d+ 2 (q~ 1 ) -go- sw -1 )? 2 


G’iq- 1 ) = 

(G Nl -d-i(q *) 90 9iQ 1 ••• 9N!-d-iq (Nl d l) )q Nl d 



( G N 2 -d-i(q x ) go giq 1 • • • 9N 2 -d-iq {N2 d 1] )g N2 d J 
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and 


(■ F d+1 (q-') 

F d+ 2(q~ l ) 

F{q ~ 1] = Fd+N^q- 1 ) 

. Fd+M^q' 1 ) 

Thus, the predicted outputs over the prediction horizon are given by 

y(* + Wi) 
y(k 4- Ni 4 1) 
y = : 

y(k 4 N 2 ) 

where, 

dNi-d-i gNx-d-2 **• t) 

dNx-d 9Ni-d-l **• 0 

9N 2 -d-l 9N 2 -d - 2 * • • 9N 2 -N u ~d 

and 

' Fd+M^q- 1 ) 4 
f = : y(k) 

. Fd+Niiq -1 ) . 

(G^-d-iiq- 1 ) -g 0 - 9\q~ l 9N l -d-\q~ {Nl ~ d ~ l) )q Nl ~ d 

4 \ A u(k — 1) 

_ (GAr 2 _d_i(? _1 ) - go — 9\9~ l 9N 2 -d-iq~ {N2 ~ d ~ 1) )q N2 ~ d . 

A.l Control law 

Using Eq. (A. 8), the cost function (Eq.(3)) can be written as 

J(Ni,N 2 ,N u , A) — (GAu 4 f- w) t (GAu 4 f— w) 4 AAu r Au (A.9) 

where 

w(t 4 d 4 N\) 
w(t 4 d 4 Ni 4 1) 

w = 

w(t 4 d 4 N 2 ) 

Equation (A.9) can be re-written as 

J(Ni, N 2 , N u , A) = iAu T HAu + bAu + f 0 (A.10) 
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where 


H = 2(G r G + AI) 
b = 2(f-w) r G 
f 0 = (f-w) T (f-w) 

Minimizing cost function (3) to obtain optimal control is an unconstrained optimization 
problem as no other constraints are considered on control input or plant output. The 
optimal input Au* is obtained by setting the gradient of J with respect to Au equal to 
zero. That is, 

Au* = — H -1 b T - (G t G + AI) -1 G r (w - f) (A.ll) 

The actual control signal sent to the plant is 

Au*(k) = JiAu* (A. 12) 

where, I\ — [1, 0, 0, ... ,0] is the 1 x N u row vector. Prom Eq. (A. 12) it can be seen that 
A u*(k) is the first element of vector Au*. 


33 



B System transformation 

Typically, most of the system equations are available in the form with the input being u 
and not A u. However, the GPC control law developed in Eq. (26) is based on the system 
having input A u. This section is intended to provide the equations to transform the system 
with input u to a system with input A u. 

Consider a MIMO LTI system with p inputs and q outputs. 

x(k + 1) = Ax(fc) -f Bu(fc) 

y (k) = Cx(k) + Du(fc) (B.13) 

where input to the system is u and output is y. Let H(z) denote the transfer function from 
u (k) to y(fc). Then, 


Y(z) = H(z)U(z) 


where H(z) has a realization 


H = 


‘ A 

B ' 

C 

D 


Now let T(z) be the transfer function from new input Au(fc) to u (fc), i. e., 

U(z) = T(z)AU(z) 

where, T(z) is p x p transfer matrix and has the form: 


T(z) = 

Suppose T (z) has a reahzation 

T = 

and the state 


l-z~ l 


l 


pxp 


' A a 

B a 


Ipxp 

Ipxp 

Ca 

D a 


Ipxp 

Ipxp 


(B.14) 


(B.15) 


xa(A:) = u (k - 1). 

Now, let II a (z) denotes the system with input Au (k), output y(fc), and having following 
state-space description. 


x(k + 1) = Ax(fc) + BAu(A:) 
y (k) = C x(k) + DAu(fc) 


Then, 


Y(z) = H(z)T(z)AU(z) 
= H a (z)AU(A:) 
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where Ha(z) has realization 


H a 


A B 
C D 


A BC a 
0 A a 
C DCa 


BD a 

B a 

DD a 


and the new state is 


x(fc) = 


x(fc) 


9.(k) 

X A (fc) 


n(k — 1) 
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C Proof of Lemma 1 


We will use the principle of optimality to derive the optimal control law. The principle of 
optimality is well known and is stated below The Principle of Optimality: An optimal 
policy has the property that whatever the initial state and the initial decision are , the re- 
maining decisions must constitute an optimal policy with regard to the state resulting from 
the first decision. 

In the context of our problem, the principle may be stated as follows: 

If a closed-loop control u*(k ) — f(x(k)) is optimal over the interval 0 < k < N, it is 
also optimal over any subinterval m < k < N , where 0 < m < N. 

Let us define a set of various costs Fj as follows 

F n = x T (k + N)(Q + C T C)x(k + N) 

Fj = x T (k + j)C T Cx(k + j) + \Au T {k + j)Au(k + j) (j = 0, 1, • • • ,N - 1) 
Then, we an re-write the cost function, Eq. (22), as follows 

Jn = Fo 4- Fi 4 b Fjv-i + F N ~ x r (A;)C r Cx(k) 

Note that x(/c) does not affect the optimal control law Au*(j), for j > k. 

Let S m (m = 1, 2, • • ♦ , N + 1) be the cost from j — k + N — m-\-l to j = k + N. Then 
S m can be written as 


S m = Jn - Jn-tu — Fjsf-m+i H ^ F n 

From the Principle of Optimality if Jn is the optimal cost, then so is S m . By the definition 
of S m , we know that 


Sjn+l — S m + Fj\f — m 

— S m + x t (A; + N — m)C T Cx(k + N — m) 4- 

AAu t (A: + iV- m) Au(k + N - m) (C.17) 


If is the optimal cost then we can get the optimal control Au*(£ 4- iV - m) by replacing 
S m by in Eq. (C.17) and then setting 


dS m +i 

<9Au (k + N — m) 


(C.18) 


We will obtain the expression for the optimal control by the method of induction. Consider 
S\. By the assumptions of GPC, Au (k + j) = 0 for j > N. So 


where P is defined as 


s* = x T (k + N)(Q + C T C)x(k + N) 
x T (k.+ N) P(k + N)x(k + N) 


P(k + N) = Q + C T C 


(C.19) 
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Prom Eq. (C.17) ; we get 
S 2 = Si + Fn - x 

= SJ + x T (fc + N- l)C r Cx(fc + JV - 1) + AAu T (k + N- l)Au(fc + N- 1) 

Substituting for and using state-space equation x(A; + 1) = Ax(fc) + BAtt(fc), we get 

S 2 = [Ax(k + N-l) + BAu(k + N-l)] T P(k + N){Ax{k + N-l) + BAu(k + N-l)] 
+x T (fc + N- 1)C T C x(Jfc + JV - 1) + AAu r (fc + JV - l)Au(fc + JV - 1) 

Now setting 

^ .0 

dAu(k + N- 1) 

we get optimal control law 

Au *{k + N- 1) = -K(k + N- l)x(Jfc 4- N - 1) (C.20) 

where 

K {k + N - 1) = [B t (Q + C r C)B + AI]- 1 B t (Q + C T C)A 

and 

SJ = x r (fc + N - 1) x 

[(A - BK(k + N- l)) r (Q + C r C)( A - BK(k + N- 1)) + 

C T C + XK T (k + N- l)K(k + N- 1)] x 
x(fc + A - 1) 

:= x T (fc JV — 1).P(& + N — l)x(k + N — 1) (C.21) 

P{k 4- N — 1) is thus given by 

P(* + JV-1) = [(A-BF(k + N-l)) T (Q + C T C)(A-BK(k + N-l)) + 

C T C + XK T (k + N — l)K(k + JV — 1)] (C.22) 

Prom Eqs. (C.19) and (C.21), by induction we can write 

= x. T (k + l)P(fe + N — m 4- l)x(fc 4- iV — m 4- 1) 

and 

Sm+i = 4- x T (fc 4- JV — m)C T Cx(k 4- N — m) 4- AAu T (fc 4- N — m)Au(/c 4JV-m) 

= [Ax(fc + N — m) 4- BAu(fc 4- N - m)] T P(k 4- JV — m 4- 1) x 
[Ax(fc 4- N - m) 4- B Au(A; 4- JV - m)} 

4 -x r (fc 4- JV - m)C r Cx(fc 4- JV — m) 4- AAu r (fc 4- iV — m) Au(fc 4- iV — m) 

Now, Eq. (C.18) can be used to find the optimal control law as follows 

= 0 =► An*(k + N-m) = -K(k + N - m)x.(h + N - m) (C.23) 

dAu(k + N — m) 
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where 


K{k + N - to) = [B T P(k + N-m + 1)B + AI]" 1 B :r P(A: + N - m + 1)A (C.24) 

Substituting for Au yields 

5 m+i = xT ( k + N -m)x 

[(A - B K(k + N- m)) T P{k + N-m+ 1)(A - B K(k + N-m)) + 
C T C + XK T (k + N - m)K(k + N — to)] x 
x(/c + N — m) 

:= x T (/c + N - m)P(k + N — m)x(/c + iV — ra) 
where, P(/c + fV — m) is given by 

P(A; 4- JV - m) = A T P(k + N-m+ 1)(A - BAX* + A - m)) + C T G (C.25) 

and the optimal cost is given by 

S*N + 1 = x T (k)P(k)yi(k) 

Jn = S* N+1 - x T (fc)C r Cx(jfc) 

= x T (k)(P(k) — C T C)x.(k) 

From Eqs. (C.20), (C.24) and (C.25), the optimal control law Au*(fc) is given by 

Au *(fc) - -A(fc)x(fc) (C.26) 

where, A(fc) is obtained from the following recursive equation 

K(k + j) = [B T P(k + j + 1)B + AI^B^A; + j + 1) A (C.27) 

and P(k) is the solution of the Riccati Difference Equation (RDE) 

P(k+j) = A T P(ifc+j+l)A-A T P(A:+j+l)B[B T P(fc+i+l)B+AI] _1 B r P(fc+j+l)A+C T C 

(C.28) 


with the boundary condition 


P{k + N) = Q + C T C 


(C.29) 
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Proof of Lemma 2 


In this section, we will show that the control law in Eq. (C.26) is same as the one in Eq. 
(26). We will show by manipulating these equations that their right hand side (RHS) and 
left hand side (LHS) are same. 

Pre-multiplying both sides of Eq. (26) by (g t G + AI + C T QC^j we get 
(G r G + AI + C r QC) Au* = - (g t L + C r QA") x(fc) 


Now, substituting G, C, and L, in above equation yields 
B r C r CB + --- + B r (A JV - 1 ) T (Q + C T C)(A iV - 1 ) r B ••• B T (A JV - 1 ) r (Q + C T C)B 


+AIAu* = - 


B t (Q + C r C)A ;v - 1 B • • - B r (Q + C t C)B 

( B t C t CA + B T A r C T CA 2 + • • - + B r (A JV-1 ) r (Q + C r C)A JV \ 
B r C T CA 2 + • - ■ + B t (A n ~ 1 ) t (Q + C t C)A n 


Au* 


V 


B r (Q + C t C)A'" 


x(A;)(D.30) 


J 


The i th equation is given by 


[b t C t CA’ - 1 B + B r A r C T CA 4 B + • - • + B r {A N ~ i ) T {Q + C 7 'C)A jv ~ 1 b] Au *(J fc) 

+ [B T C r CA i - 2 B + B T A r C r CA* _1 B + - - • + B r (A^ _i ) r (Q 4- C r C)A yv_2 B] Au*(fc + 1) 


+ 

+ 

+ 

+ 


B t C t CB + B r A r C r CAB + • • - + B r (A ;v - i ) T (Q + C r C)A iV - < B] Au*(jfe + i - 1) 
B r A T C T CB + B T (A 2 ) T C T CAB + - • • + B r (A N - < ) T (Q + C t C)A jv - < - 1 b] Au*(fc + i) 


+ [B T (A ;v - i ) T (Q + C T C)B] Au *{k + N- 1) 

-[* 


B r C r CA* + B 2 A 2 C 2 CA t+1 + • • • + B 


T A Tr'Tf' A *+l 


Now from Eq. (C.26), we can write 


r (A N ~ i ) T (Q + C r C)A JV ] x(Jb) (D.31) 


Au*(Hi) = —K (k + i)x(fc + i) 

= - [B r P(& + i + 1)B + Al] _1 B T P(fc + i + 1)A x(* + i) 

Pre-multiplying above equation by j^B T P(fc + i + 1)B + Alj gives 

[B r P(fc + i + 1)B + Al] Au*(fc + i) = -B T P(k + i + l)Ax(A; + i) 

. . . (* = 0, 1,2, • • • , JV — 1) (D.32) 

Now, using 


x(fc + i) 


i- 1 

= A l x(fc) -f ^2 A t ” J "" 1 BAu(A: + j) 


o 


(D.33) 
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= A i x(fc) + [ A i-1 B A i_2 B • • • B ] 
we can write Eq. (D.32) as 

B r P(k + i + 1) [ A J B A' _I B • • • A‘ _1 B B ] 


{ Au(/c) \ 
Au(fc 4* 1) 

Au(£; 4- i — 1) J 


(D.34) 


( Au *(*) \ 

Au*(k + 1) 

^ Au *(k + i) J 


+XIAu*(k + i) = -B T P{k 4 i + l)A i+1 x(fc) (i = 0, 1, 2, • • • , N - 1) 

Substituting for index i , we get 

\ ( Au*(fc) \ 
Au*(fc 4 1) 


( B T P(k 4 1)B 

B T P{k 4 2) AB B T P(k 4 2)B 


V B T P(k 4- N)A n ~ 1 B B T P{k + N) A N ~ 2 B ■ ■ - B T P{k 4 N ) B J \ Au *(k 4- N - 1) / 


AU* 


+AIAu* = - 


( B T P(k 4 1)A \ 
B T P(k 4 2)A 2 

V B T P(k 4 N)A n J 


x(fc) 


(D.35) 


Note that Eqs. (D.30) and (D.35) have similar structure. Now, if we transform the RHS 
of Eq. (D.35) to match that of Eq. (D.30) and, as a result, if we get the LHS of both 
equations to match we have proved the lemma. 

So, first transform the RHS of Eq. (D.35) to match that of Eq. (D.30). Consider the 
i th row in the RHS of Eq. (D.35). By Eq. (C.25), we have 

P{k + i) = A T P(k + i + l)A + C T C-A T P(k + i + l)BK(k + i) 

— A T [A T P(k 4 i 4 2)A + C^C — A ^ P{k 4 i 4 2)BA(fc 4* i 4 1)]A 4 C^C 
-A T P(k 4 i 4 1)B K(k + i) 


= (A N ~ i ) T P(k 4 N)A n ~* 4- (A^ -i-1 ) :r C' r CA JV_i-1 + (A N ~ i - 2 ) T C T CA N - i - 2 4 
■ • • 4 C T C - (A N - i ) T P{k 4 N)BK(k + N- 1) A*-* -1 
—(A N ~ i ~ 1 ) T P(k + N - 1)B K(k + N- 2)A N ~ i ~ 2 

-A T P(k 4 i 4 1)B K{k 4 i) (D .36) 

That is, the i th row of RHS in Eq. (D.35) becomes 

-B r P(fc 4 i) A* - -$ + D 

where 

$ = [B r (A JV_i ) r P(A: 4- N)A n 4- B T (A Ar ~ i-1 )' r C T CA Ar_1 4- ¥ B T C' r C.A(ijk$£) 
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and fi = \B T (A N ~ i ) T P(k + N)BK(k + N - l)A N -‘ (D.38) 

B T (A N ~ i ~ 1 ) T P(k + N- 1)B K{k + N- 2)A N ~ 2 (D.39) 

B T A T P{k + i + 1)B K(k + i)A i ]x(fc) (D.40) 

Note that, # is same as the RHS of Eq. (D.31). Next, we will move Q to the LHS, and 
then simplify the equation. Now, by (D.33), we have 

i — 1 

A'xflfe) = x(fc + i) - A i_J_1 BAu(A: + j) (D.41) 

j — 0 

(D.42) 

Substituting Eq. (D.41) in Q and using optimal control law, Au*(k+j) = —K(k+j)it(k+j), 
we get 

n = 

—B T (A N ~ i ) T P(k + JV)BAu *(fc + N- 1) 

-B r (A x-tfPik + N)BK(k + N- 1) (a n_2 BAu (k) + • • • + BAu(fc + N- 2)) 
—B T (A N ~ i ~ 1 ) T P(k + N- l)BAu *(k + N- 2) 

-B r (A- v ~ i_1 ) T P(A; + N- 1 )BK{k + N- 2) (a jV_3 BAu(A;) + • • • + BAu(k + N - 3)) 
-B r A r P(Jb + i + l)BAu*(fc + i) 

-B T A T P(k + i + l)BK(k + i ) (A i-1 BAu(fc) + • • • + BAu(A + i - 1)) 

The the ith row of LHS of Eq. (C.26) is 

B T P(k + i)A i ~ 1 BAu*(k) + B r P(fc + i)A i_2 BAu*(Jfc + 1) • • • 

+ B T P(k + i)BAu*(fc -M — 1) (D.43) 

Substituting Eq. (D.43) in Eq. (D.36) and moving ft to the LHS, we get the coefficients 
of Au*(k) on LHS as 

B T [(A 7V - i ) T P(fc + AT)A n-< + (A N - i - 1 ) T C T CA Ar - i - 1 + • • • + C r C]A i - 1 B 

where, P(k + N) = C T C + Q. This coefficient is same as the one in Eq. (D.31). Similarly, 
we can relate the coefficients of Au*(fc + j) (J = 1, 2, • • • , i — 1) to the ones in Eq. (D.31). 
Next, look at the coefficients of Au *(k + j), (j = i, i + 1, • • • , N — 1). For j = N — 1, the 
coefficient of Au ’(k+j)) on LHS becomes B T (A N ~ t ) T P(k + N)B, which is same as in Eq. 
(D.31). For j = i, the coefficient of An*{k + j)) on LHS is given by 

B t A T P(k + i + 1)B + B T (A 2 ) T P{k + i + 2)B K(k + i + 1)B 
+ • • • + B r ( A^ - * -1 ) r P(k + N — 1)B K(k + N — 2)A^- 3 B 
+B T (A jV - < ) T P(fc + N)BK(k + N- 1)A JV - <_2 B 

From Eq. (D.36), we have 
P(k + i + 1) 

= (A^-'f^/t + JVJA^- 1 ^ (A Ar - i - 2 > : T.cICA- Jv rL-A+, • • _+ cTc_ 

— (A JV-i-1 ) r P(fc + N)BK{k + N- 1)A n ~'~ 2 

—(A N ~ i ~ 2 ) T P(k + N- 1)B K{k + N — 2)A N ~ i ~ 3 

-A T P(k + i + 2)B K(k + i + 1) 
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Substituting P(k 4- i -f 1), we get the coefficient of Au *(k + i) as follows 
B T (A N ~ i ) T P(k + N)A n ^ x B + B T (A jV - i - 1 ) T C T CA Ar ’ i “ 2 B 

+ • +b t a t c t cb 

which is same as the one in Eq. (D.31). 

Thus, we have shown that the optimal control law given by Eq. (26) is same as the one 
given by Eq. (C.26) ■ 

Numerical Example: 


Consider an LTI system, 


where 


x(fc + 1) = Ax(fc) + BAu (k) 
y(k) = Cx(fc) 


A = 


1 2 
3 1 




and the cost function of Eq. (22). Let the GPC control parameters be A = 1, N\ = 1, 
N 2 — N u — 3, and 


Q = [!] • 


Now the GPC control law computed by Eq. (26) yields 


K t 


gpc 


(G t G + AI + C r QC)- 1 (CQA 3 + G t L ) 

' 3.6368 2.4684 

1.8240 3.1017 

0.1072 -0.3723 


Similarly, the LQ control law computation for the same cost function using Eq. (D.35) 
yields 


Klqr 


B t P(1)B + A 

B r P(2)AB B r P(2)B + A 
B t P(3)A 2 B B t P(3)AB B t P(3)B + A 


B r P(l)A \ 
B t P(2)A 2 
B t P(3)A 3 ) 


3.6368 2.4684 

1.8240 3.1017 

0.1072 -0.3723 


which is same as K gpc . This shows that both control laws are the same at time instant 
t = k. When the cost horizons for GPC and LQ strategies are identical. 
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